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Dynamics of one-dimensional systems of Brownian particles with short-range repulsive interaction 
(diameter a) is studied with a liquid-theoretical approach. The mean square displacement, the 
two-particle displacement correlation, and the overlap-density-based generalized susceptibility are 
calculated analytically by way of the Lagrangian correlation of the interparticulate space, instead 
of the Eulerian correlation of density that is commonly used in standard mode-coupling theory. In 
regard to the mean square displacement, the linear analysis reproduces the established result on the 
asymptotic subdiffusive behavior of the system. A finite-time correction is given by a Lagrangian 
version of mode-coupling theory. The notorious difficulty in derivation of mode-coupling theory 
concerning violation of the fluctuation-dissipation theorem is found to disappear by virtue of the 
Lagrangian description. The Lagrangian description also facilitates analytical calculation of four- 
point correlations in the space-time, such as the two-particle displacement correlation. The two- 
particle displacement correlation, which is asymptotically self-similar in the space-time, illustrates 
how the cage effect confines each particle within a short radius on one hand and creates collective 
motion of numerous particles on the other hand. As the time elapses, the correlation length grows 
unlimitedly, and the generalized susceptibility based on the overlap density converges to a finite value 
which is an increasing function of the density. The distribution function behind these dynamical four- 
point correlations and its extension to three-dimensional cases, respecting the tensorial character of 
the two-particle displacement correlation, are also discussed. 



I. INTRODUCTION 



In spite of the ever more active studies on glassy sys- 
tems in the last few decades [l|, 0, and even after the 
discovery of dynamical heterogeneity 0, H|, theoretical 
understanding of jamming, plasticity, and glass transi- 
tion remains a challenging problem in statistical physics. 
The dynamical heterogeneity in glassy systems refers to 
the presence of non-trivial collective motion of numer- 
ous particles with some correlation length and lifetime; 
the growth of this dynamical correlation length, as op- 
posed to the static correlation which has been sought 
fruitlessly, is now considered to account for the slow- 
down of the dynamics near the glass transition. However, 
kinetic-theoretical treatment of this collective motion is 
formidably difficult, as its correct description requires a 
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four-point correlation function, whose analytical calcula- 
tion has been infeasible except for some special cases such 
as kinetically constrained models on a lattice Indeed, 
a growing dynamical length in supercooled liquids was re- 
cently calculated [H] by extending a prevalent theoretical 
scheme known as the mode-coupling theory (MCT) [3], 
but MCT suffers from its own difficulties Besides, 
it is not clear how the cage effect, which confines each 
liquid particle within a short radius, is related to the 
long-ranged dynamical correlation. 

In such a state of affairs, what Hilbcrt said in his well- 
known speech [Io| is thought-provoking. Emphasizing 
the vital role of good problems that serve to temper the 
sword of mathematical methodology, he pointed out the 
importance of proper specialization: in order to overcome 
mathematical difficulty, we should study "problems sim- 
pler and easier than the presented one" in depth. In 
solving the specialized, "simpler and easier" problems, 
one has a chance of developing an innovative method 
with gcncralizable concepts, which eventually makes it 



2 



possible to overcome the original difficulties. 

In our case, the "presented problem" is to elucidate the 
nature of the slow dynamics characterizing glassy sys- 
tems, by extracting an effective theory for it from their 
particle-scale models. To be specific, let us consider dense 
colloid suspensions, modeled as a three-dimensional sys- 
tem of interacting Brownian particles: 

d 



A t *i-3^X^fo*)+A*(*) 5 (1-1) 



j<k 



(1.2) 



Without the interaction (V = 0), the statistically aver- 
aged or coarse-grained density field would be governed 
by the simple diffusion equation, with the diffusion con- 
stant D = k^T/p. Contrastivcly, for repulsive short- 
range interaction potentials with excluded volume effect 
that determines the particle diameter a (and with some 
contrivance to avoid crystallization), the Monte Carlo 
and molecular-d yna mical equivalents of Eq. (jl.ip show 
glassy behavior [ill \v\ . As the density increases, the 
cages have stronger effect, which results in the extremely 
slow relaxation of the intermediate scattering function 
F(q, t) oc (p(q, t)p(— q, 0)), with p denoting the Fourier 
component of the density, such that 



-iq r 



(1.3) 



The slowdown of the dynamics suggests the presence of 
some spatial structure with growing scale, but this is not 
detectable via the static structure factor, S(q) — F(q, 0). 
Instead, a cooperative motion with growing correlation 
length is evidenced by four-point dynamical correlations 

Hani, such as 



X4(t) 



(p(r 1 ,t)p(r 1 ,0)p(r 2 ,t)p(r 2 ,0)} 



-<p(n,t)p(n,0)) (p(r a ,t)p(r 2 ,0)) 



d 3 rid 3 r 2 . 



In particle-based computations, it is evaluated typically 
via the overlap density 



as 



X 4(t) « <Q(*) 2 ) - m))\ 



(1.4) 



with S a denoting the overlapping function, which has a fi- 
nite value around r = and vanishes for r»o. The spa- 
tiotemporal structure, characterized by four-point space- 
time correlations, is referred to as dynamical heterogene- 
ity and now studied intensively Q, in the hope of find- 
ing the secret of the glassy dynamics in it. Neverthe- 
less, calculation of dynamical four-point correlations for 
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FIG. 1: The mean square displacement (R 2 ) in SFD, plotted 
as a function of t. The dotted line indicates the slope for 
i 1/2 . The system size is specified as N = 2 15 = 32768 and 
po = N/L = 0.25 a~ 1 . In the inset, the same (-R 2 ) is plotted 
against \ft for < \[Dt < 10 a. 



Eq. (jl.ip is so difficult that, to our knowledge, no one 
has managed to calculate them analyticall y. The theo- 
retical calculation of \4 by Donati et al. [la ] concerns 
the static case. The most remarkable success was made 
by Toninclli et al. Q, who calculated Xi analytically 
for elastic solids and kinetically constrained lattice mod- 
els. These are cases in which permanent "bonds" are 
present: the situation is different in the colloidal system 
governed by Eq. (jl.ip , which lacks such "bonds" connect- 
ing the particles (at least explicitly). Besides, it should 
be emphasized that MCT 0, [|[ targets on the dynamical 
structure factor F(q, t) and not on four-point correla- 
tions. As long as the standard variables such as p(q, t) 
in Eq. (| 1 .3[) are used, a four-point correlation function 
implies a four-body correlation. Since MCT is a closure 
theory in which quadruple (four-body) correlations are 
approximated by products of F, it is unlikely to describe 
four-point correlations accurately, for example, with the 
accuracy needed for treatment of entanglement effects in 
melts of rod polymers [TB]. Although a calculation of 
three-point correlation within the MCT approximation 
was recently reported still X4 remains insurmount- 
able. 

With this situation in mind, according to Hilbert's ad- 
vice, let us look for "simpler and easier" problems that 
may give insight into real glassy systems. Instead of 
Eq. (jl.ip . here we study its one-dimensional version: 

mXi = -pXi - ^7 E v[Xk ~ X ^ + (1 ' 5) 



j<k 



whose behavior is known by the name of single-file dif- 
fusion (SFD) [1714251 ]. For V we adopt a short-ranged 
repulsive potential, such as Eq. (|A1|) in Appendix [A] 
so that the system exhibits a liquid-like structure factor 
[H, H3 ■ Despite this retainment of the original feature 
of Eq. (jl.ip . since the one-dimensional system does not 
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exhibit glass transition, Eq. (|1.5j) may seem less than a 
toy model for studies of glassy dynamics. This is a su- 
perficial view, however, if we note that SFD throws one 
of the difficulties of MCT into relief — the difficulty met 
by Miyazaki and Yethiraj [l6| in their study of rod poly- 
mers, for which they needed to go beyond the standard 
MCT by directly calculating the four-point correlation, 
as will be explained later. 

It is also a kind of four-point correlation function that 
we treat here to account for the slow diffusion subject 
(11.5 1. By "slow" we mean that SFD is sub- 
19 . l2(il l: by observing longtime behavior of the 



to Eq 
diffusive 

mean square displacement (MSD), denoted by (R 2 ) with 



Rj = Rift) = X 3 (t) - Xj(0), 



(1.6) 



it is found that (Rj) for SFD behaves like \ft [see 
Eq. (|1.7p shown later], which is slower than (Rj) oc t 
expected for the normal diffusion. The subscript in (Rj) 
can be omitted as we assume that the system is statisti- 
cally uniform. The subdiffusive law, (R 2 ) oc ^/i, is read- 
ily confirmed by simulation of a system with TV particles 
in a periodic box of the size L. In the calculations shown 
in Fig. [TJ the mean density, po = N/L, equals 0.25 er -1 , 
and the system is statistically steady. For computational 
details, see Appendix [3] The finite size effect is elimi- 
nated by taking sufficiently large L and interpreting the 
word "longtime" as p^ 2 <C Dt <C L 2 —> 00. For the 
longtime regime in this sense, Kollmann showed [20I ] 



/ i? 2) = ^./^ cx< l/2 j 
PO V 7T 



(1.7) 



where S = S(0) denotes the long wave limiting value of 
the static structure factor S(k), and D c is the collective 
diffusion constant. We note that Kollmann also needed 
to calculate a kind of four-point correlation (k' 2 ' in his 
notation) in derivation of Eq. (| 1 . T[) . 

The slowness of SFD is ultimately due to the presence 
of the repulsive potential term in Eq. (|1.5[) . The structure 
of Eq. (|1.5j) is basically the same as that of the three- 
dimensional equation (jl.l[) . except for the geometrical 
difference that makes the cages temporary in the three- 
dimensional systems while the one-dimensional cages are 
eternal. Then there is a question: Docs a straightfor- 
ward application of MCT to Eq. (|1.5p reproduce the sub- 
diffusive law (|1.7[l ? Unfortunately, the answer is neg- 
ative according to Miyazaki [§, [28|, who demonstrated 
that MCT predicts normal diffusion, (R 2 ) oc t, instead 
of t 1 / 2 . His results imply that the dynamical constraint 
by the long-lived cages is not described accurately enough 
by the standard MCT. Essentially the same difficulty oc- 
curs in MCT for rod polymers (l6| . as it predicts that 
the self-diffusion coefficient becomes isotropic for large 
aspect ratio and therefore fails to describe the entangle- 
ment effects. 

Motivated by Miyazaki's study on SFD [HUll and also 
inspired by the notion of Lagrangian correlation in clo- 



sure theories of fluid turbulence [29h32| . the present au- 
thors reformulated the continuum theory of SFD with 
the label variable (24|. This is essentially an application 
of the Lagrangian description [H, HH , as opposed to the 
Eulerian, to the Langevin equation for the density field. 
The label variable method made it possible to reproduce 
the subdiffusive law (|1.7[) . with the correct power and 
correct coefficient. 

In the present paper, we demonstrate that the label 
variable method allows us to calculate four-point corre- 
lations explicitly It is demonstrated by generalizing the 
calculation of MSD i24j to the two-particle displacement 
correlation (2pDC), 

(RiRj) = ([Xi(t) - Xi(p)] [Xj(t) - Xj(0)) ) . (1.8) 

Quantities analogous t o 2p DC have been studied by a 
number of researchers |l2l . I35T437T ] with numerical data 
from molecular dynamics, and it is also the main ingre- 
dient of the theory of \i fo r elastic waves by Toninclli 
et al. Q. Nevertheless, analytical calculations of 2pDC 
for "bondless" particle systems have never been reported. 
Here we calculate 2pDC analytically in terms of gener- 
alizablc and liquid-oriented concepts, so that the theory 
could be extended to truly bondless systems in the near 
future. At first, we calculate 2pDC as a function of the 
elapsed time t and some properly defined label distance 
which coincides with \i — j\ in the one-dimensional cases; 
and subsequently, we show it to be re-expressible as a 
function of t and the initial distance d — Xj(0) — Xi(0), 
which we denote with Xn(d,t). For i = j, Eq. (|1.8|) re- 
duces to MSD. The two-particle correlation, (RiRj) with 
i j, provides an intuitive form of four-point correlation 
in comparison to Xiy the 2pDC with i =/= j vanishes for 
free Brownian particles and, for SFD, evidences the clus- 
ter size that behaves like po y/D c t, accounting for the slow 
diffusion. In addition, from the knowledge of (RiRj), we 
can perform a fully analytical calculation for the self part 
of X4, defined as 

x!W«(Qi)-(Qs) 2 , 

where Q s = Eja(Ri(t)) [see Eq. in Sec. [V] for 

details]. Reflecting the eternity of the one-dimensional 
cages, the longtime limiting value, X4(+°°)> is finite, as 
will be shown in Eq. (|5.17|) . 

Furthermore, we demonstrate that an MCT-like non- 
linear theory for the fluctuation of 1/p can be devel- 
oped in a field-theoretical style, without violating the 
fluctuation-dissipation theorem (FDT). This is possible 
because the problem of the multiplicative noise is natu- 
rally resolved by virtue of the label variable. The contri- 
bution of the memory term gives a finite-time correction 
to Kollmann's asymptotic result (|1.7[) , visible as a finite 
intercept of the asymprotic straight line on the v^-axis 
in the inset in Fig. [TJ 

The paper is organized as follows: After summariz- 
ing in Sec. [II] the idea of the continuous label variable 
method and some of its results, we apply it to the cal- 
culation of the two-particle displacement correlation in 



Sec. Mil Strictly speaking, what we present in Sec. Mil is 
not (RiRj) itself but its continuum equivalent, calculated 
theoretically for the longtime regime. Subsequently, in 
Sec. lIVI we demonstrate a systematic and FDT-preserving 
derivation of an MCT-like equation in the Lagrangian 
description. The "Lagrangian" MCT equation provides 
us with the finite-time correction to MSD and the two- 
particle displacement correlation. With this finite-time 
correction taken into account, two different forms of four- 
point correlation functions are calculated in Sec. [Vj one 
is Xft(d,t) and the other is xf (*)• Section [VI] is allotted 
for discussion and concluding remarks. 



II. CONTINUUM THEORY OF SINGLE-FILE 
DIFFUSION 

Let us begin with summarizing our previous results on 
MSD for SFD J24|. By "continuum theory" we mean that 
the theory is formulated in terms of some hydrodynamic 
quantity such as the density p, rather than direct treat- 
ment of the particles. Our idea consists in adoption of 
the continuous label variable £, which we take instead of 
the position x as the independent variable, and we also 
change the dependent variable from the density p(x, t) to 
the fluctuation of the particle interval, denoting it with 
tp(£,t). On the basis of the correlation of tp calculated 
for the longtime regime, we can re-derive Eq. (| L . T[) . 

Since our label variable method is intended as a re- 
formulation of MCT, we start from essentially the same 
Langevin equation as in field-theoretical formulation of 
MCT for dense colloidal suspension 38[. The Langevin 
equation, derived from Eq. (|1.5[) for the density p(x, t) = 

YjjPj witn Pj =S(x- 
is given as follows: 



Xj(t)) and its flux Q = y\ PjXj 



d t p + d x Q = 0, (2.1a) 
Q = -D (d xP + T^d x U] +5> i (s > *)/,•(*), (2.1b) 



U = U[p](x)= dx'V{x-x')p{x'). 



(2.1c) 



By eliminating Q and introducing 

f p {x,t) = -d x ^2pj(x,t)fj{t), 

3 

we write down the equation for p(x, t) as 



d t p(x, t) = Dd x \d xP + d x Uj + f p (x, t) (2.2) 

with d x ll = d x U[p](x) and 

(f p (x,t)f p (x',t')) = 2Dd x d x >p(x,t)5(x-x')5(t-t'). 

(2.3) 

Note the presence of p(x, t) on the rig ht-hand side of 
Eq. (|2.3|) : the noise is multiplicative [33 |. 




FIG. 2: (Color online) A schematic view of the mapping from 
£ to x, given by inverting the function £ = f (x, t) in Eq. (|2.6I) . 
with t fixed arbitrarily. The label and the position of the ji-th 
particle is denoted by Sj and Xj, respectively. 



As the particles have excluded volume effect and there- 
fore cannot overlap, the barrier expressed by V must be 
infinitely high. In the one-dimensional system, this bar- 
rier acts as a topological constraint, or the "no overtak- 
ing" rule in plain words, which keeps the order of the 
particles. In MCT, however, V is replaced with a fi- 
nite effective potential, either explicitly by hand (in the 
diagrammatic formalism) or implicitly by the projection 
operator (in the Mori-Zwanzig formalism). As a result, 
MCT cannot incorporate the "no overtaking" effect of 
V properly. An asymptotic analysis of one-dimensional 
MCT for longtime limitgives the normal diffusion, con- 
tradicting to Eq. (JTTTJ) [a Hi]. 

When Miyazaki gave a talk on SFD (28|, the present 
authors were nursing an idea of importing the concept 
of Lagrangian correlation from the theory of fluid tur- 
bulence |3(j-[32| hrto the studies of glassy systems, and 
this idea was soon noticed to be applicable to SFD. The 
idea consists in adoption of the continuous label variable, 
which wc denote with £. 

Though it is popular in continuum mechanics to take 
the initial position of each material element to label it, 
here we define £ in a different way, avoiding to trace the 
whole history of the system back to the initial configura- 
tion. Wc construct the label function £ = £(x, t) so as to 
satisfy the following three requirements: 

1. The label should satisfy the convective equation 



(pdt + Qd x )£(x,t)=0. 



(2.4) 



2. The label should be related to the snapshot of p 
and Q in such a way that the continuity equation 
(I2.1al) is satisfied. 

3. The function needs to be invertible, in the sense 
that a mapping from (£, t) to x = x(£,t) should 
exist. 
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To satisfy the second requirement, we utilize the fact 
that Poincare's lemma [4(j is applicable to the continuity 
equation (|2.1a|) , which guarantees the existence of £ such 
that 

d x £(x,t) = p, 8tt(x,t) = -Q; (2.5) 
in fact, a solution to Eq. (|2.5[) is explicitly given by 

£ = t) = [ p{x\ t)dx' + const. (2.6) 

Jx a (t) 

Then it is straightforward to verify that all the three 
requirements are satisfied. To be precise, since p con- 
sists of a sum of delta functions, the integral in Eq. (|2.6[) 
gives a multiple step function (see Fig. [5]), which needs to 
be slightly smoothed to justify the single- valuedness and 
the invertibility. Although this smoothing could be in- 
terpreted physically as a consequence of coarse-graining, 
here we refrain from involving ourselves with such a deli- 
cate matter and regard the smoothing simply as a math- 
ematical rcgularization. We define 3 3 - = £(Xj(t),t) for 
j € Z, taking it for granted that the particles are num- 
bered consecutively; then, with the smoothing and the 
integral constant tuned appropriately, we have Sj = j. 

Using the label variable £ as the spatial coordinate, 
now we rewrite the Langevin eq. (|2.1[) . The chain rule 
for the differential operators gives 



d, 



dx 



P d s , (d t ■ ) x = (d t ■ ) e -Qd e . 



As a kinematic relation that replaces Eq. (|2.1a|) , from the 
identity dtd^x = d^dtx we find 



P&t) 



(?) 



(2.7) 



Then, rewriting Eq. (|2.1b[) with and substituting it 
into Eq. (|2.7[) . we obtain 



dt 



P{t,t) 



deU 



d s j2 s (t-Zi)fj(t) 



= -Ddt [d s p + 2 sinh (pad 6 ) p] + t), 

(2.8) 

where V is replaced with the effective potential as before, 
and /l satisfies 

(h(u)h(e,t')) = 2Dd^j2 s ^-^m-e)s(t-t'). 

(2.9) 

Having changed the independent variables from (x, t) 
to (£, t), we change the dependent variable as well. In- 
troducing the fluctuation of the particle interval, 



1> = 1>fot) 



Po 



(2.10) 



and its Fourier modes 



we rewrite Eq. (|2.8[) in the form 

d^(k, t) = -D^ [\ + 2sin ^° gfc ) i>(k, t) 

+ E V^(-p,^(-9,t) + 0(^ 3 ) 

k+p+q= " +Pok(k,t) (2.12) 

with D* = PqD and 

k p q 

= D*k 1 H sinpocfc + t~ s'mpoap + - — sinpocq 

\ pq kq kp 

(2.13) 



As for the statistics of the random force term, Eq. (|2.9j) 
is re-expressed as 

9D 

pl(h(k,t)f L (-k',t')) = -^±k 2 5 kk ,8{t~t')- (2.14) 

see Endnote [4l|. We also note that the linearization of 
Eq. (|2.12j) coincides with the one-dimensional version of 
Edwards- Wilkinson equation [4^, HU ■ 

To calculate MSD, we developed a formula for it in 
terms of the correlation (ip(£,t)ip(£',0)), or its Fourier 
transform (ip(k, t)tp{— k' , 0)). In Fourier representation, 
the formula reads as follows: 

W.^fSMI^ ( , 15) 



k 2 



where 144 



A 



C(k,t) = j^(^(k,t)i>(-k,0)). (2.16) 

In Sec. HIH this formula will be re-derived as a special case 
of Eq. (|3~2|) . Wc note that Eq. (|2~T5l) may be regarded 
as a refinement of the approximate formula by Alexander 
and Pincus fl9j |. 

<#) * const, x f Ml%! d , (2 . 17) 



and therefore Eq. (|2.15|) may be referred to as the modi- 
fied Alexander-Pincus formula. In the linear (Edwards- 
Wilkinson) case, a rid-dimcnsional version of Eq. (|2.17[) 
has appeared in the literature H^, |4(|; we will dis- 
cuss it later in Sec. I VII calling attention to some delicate 
points about the extension to the n<j-dimcnsional liquid 
dynamics. 

Once the formula (|2.15[) is derived, all depends on the 
knowledge of C. In particular, the long-time behavior of 
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MSD is found by linear analysis of Eq. (|2.12[) . From the 
linearized equation, the correlation C is calculated as 



COM) 



j£_ -{D,/S)k 2 t 

L 2 



(2.18) 



with the aid of Eq. (|2~Ti|) . Substituting Eq. (|2TT8"|) into 
Eq. (|2.15[) . we can successfully re-derive the subdiffusion 
law in Eq. (|1.7[) . A refinement of Eq. (|2.18p . including 
the effect of the nonlinear term, will be discussed later in 
Sec. El 



III. FOUR-POINT CORRELATION 

A. Cooperative motion in SFD 

The slow dynamics in SFD is associated with collec- 
tive motion of the particles. This collective motion, ob- 
tained by numerical integration of Eq. (|1.5[) , is depicted 
in Fig. [3J In the numerical calculations (see Appendix [S] 
for details), after the system has reached the thermal 
equilibrium, we choose some instance as t = and record 
the "initial" position of each particle, say, Xi(0). To pro- 
duce Fig. ODa), at t = 2 n x 10 a 2 /D (with n = 1,2, . . .) 
we measured the displacement Ri(t) = Xi(t) — Xi(0) for 
each i. If Ri(t) > 5er, we mark the position of the particle 
with a circle (O); if Ri < — 5cr, we mark it with a cross 
(x). As the time difference t increases, a string of the 
same kind of symbol is formed, expressing a cluster of 
particles in cooperative motion. While Fig.^a) presents 
a close-up for a relatively limited time span, a long shot 
up to t = 10 6 a 2 /D is shown as Fig. EJJb). Formation 
of large clusters in cooperative motion is visible. By the 
time difference t = 10 6 a 2 /D, a typical cluster reaches the 
size of several hundred particles, occupying a length on 
the order of 10 3 er (as each particle is assigned a space of 
1/po = 4a). 

Quantitative description of this collective motion re- 
quires some four-point correlation functions, such as 
those used for the analysis of dynamical heterogeneity 
in glassy systems. In regard to SFD, it seems espe- 
cially natural to consider the four-point correlation to 
deal with the topological constraint of the "no passing" 
rule, as the requirement that two world lines should not 
intersect involves four points in the space-time, namely 
Xi(0), XM)), Xi(t), and Xj(t). This is why Miyazaki & 
Yethiraj [16| needed to calculate a four-point correlation 
[G(l, 2; 3, 4) in their notation] to study the entanglements 
of rod polymers within the framework of liquid state the- 
ory. 

Here we demonstrate that the label variable method 
is also capable of calculating a kind of four-point corre- 
lation, in the form of two-particle displacement correla- 
tion. Though it was originally presented as (RiRj) in 
Eq. (jl.8p . there is no trouble in replacing the particle 
numbering with the label variable. Then, our task is 
derivation of a formula to calculate (R((,,t)R(t;' ,t)) with 



R{£> *) = x (£i t)~ x (£,i 0)) which will be shown in the next 
subsection. 



B. Two-Particle Displacement Correlation 

Aiming for an analytical expression for (i?(£, i)i?(£', i)) 
in SFD, let us extend Eq. (|2.15p so that we can calcu- 
late (i?(£, i)i?(£' 7 t)) from the correlation C defined by 
Eq. (|2.16[) . We start with noticing that Q/p stands for 
the velocity, whose integral in regard to t gives the dis- 
placement of the particle labelled with £: 



' p(W) 



Into this equation we substitute Q/p = dc 1 dt(l/p), ob- 
tained from Eq. (|2.7|) upon integration over £, to find 



m,t) = dt i 



i 

P0 



l + lp 

Po 



^~ik $(M)-&M)] ■ (3.1) 



Subsequently, we multiply Eq. (|3.1[) by its duplicate 
with (£, k) changed to (£',—&')) an d take the statistical 
average. The double summation on the right-hand side 
reduces to the single one, if we assume that the contribu- 
tion from the terms with k =/= k' vanishes. This is true in 
the linear case, and also seems to be justifiable for nonlin- 
ear cases within the framework of the direct-interaction 
approximation (explained later). Thus we obtain a for- 
mula allowing us to calculate 2pDC from (iptp)'- 



(R(t,t)R(e,t)) 

L 4 ro ° 



ttN 2 



dk -^0«_«, (3 . 2 ) 



with C(k,t) defined by Eq. ([2~TTJj) . Note that Eq. 
includes the modified Alexander-Pincus formula (|2.15[) 
as a special case with £ = as it ought to be. In this 
sense, the formula Q3.2p could be referred to as the ex- 
tended Alexander-Pincus formula. We emphasize that 
Eq. ()3.2p docs not rely on smallness of deformation, nor 
it requires such kind of approximation at all, as long as 
the Lagrangian description is strictly followed. 



C. Calculation of 2pDC: linear approximation 



For C in Eq. (|2 . 18[) calculated from the linear approx- 
imation of Eq. (|2.12l) . the extended Alexander-Pincus 
formula (|3.2p gives 



25 
Po 



D c t 



exp 



4p 2 D c t 



^l^lerfcJiz^ 



Po 



(3.3) 



7 



(a) 



(b) 



700 




1650 

position x/a 



1750 




-20(7 



-20cr 



FIG. 3: (Color online) Clusters in cooperative motion visualized in space-time diagrams of particles in SFD, calculated for 
po = N/L — 0.25 a -1 , (a) Worldlines of particles in the (x,t) -plane. The symbols Q and x mark particles displaced (by 
more than 5cr) rightward and leftward, respectively. The unmoving particles are indicated with small triangles. Each cluster 
in cooperative motion, involving 5 particles at least, is highlighted in a box. (b) The displacement R depicted in gray scale 
(blue-to-red scale online) as a function of £ and t. Unmoving particles are shown in white. 
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FIG. 4: (Color online) Comparison of Eq. (|3.3I) with the 
computed 2pDC in the same system as in Fig. [1] (N = 
2 15 = 32768, po = N/L = 0.25 a -1 ). The solid (red) line 



shows the self-similar solution in Eq. (|3.3I) . while the simu- 
lation data are plotted with symbols: solid triangles (a) for 
t = W0a 2 /D, open circles (O) for t = 200 a 2 / A crosses (x) 
for t = 400 a 2 /D, open squares (□) for t = 800ct 2 /D, and 
solid squares (■) for t = 1600 a 2 /D. 



This is expressible in terms of a similarity variable 

9= S i_ 
2p y/D~t 

(R(£,t)R(e,t)) 2S ( e 



(3.4) 



a.s 



- |0| erfc | 



(3.5) 

From this similarity solution we can read the size A of a 
cluster in a cooperative motion, 



A = \{t) 



(3.6) 



We have already seen such clusters in Fig. [31 though care 
should be taken in regard to the difference that Eq. (|3.6[> 
is a statistical law while Figs.[3ja) and (b) present a single 
run, and only a small portion of it is shown in Fig. G2a) . 

Equation (|3.3|) is compared with statistical results of 
particle simulations in Fig. [4] Except for the transient 
behavior slightly visible for t = 100 a 2 /D, all of the simu- 
lation results are consistent with Eq. (|3.3[) . The transient 
behavior can be studied by taking the nonlinear terms in 
Eq. (|2.12p into account, which will be discussed in the 
next section. 

The variable for the horizontal axis of Fig. [4] requires 
some consideration. As a quantity corresponding to £— 
here we have taken the distance in the label numbering, 
say, \i — j\ between the i-th and j-th particle. Though 
adequate in the present case, this is not convenient for ex- 
tensions to multidimensional cases [H, H3] , because the 
result may depend on the way of numbering. A rea- 



sonable alternative in the particle simulation is the ini- 
tial distance such as Xi(0) — Xj(0), which can be com- 
pared with the theoretical prediction by assuming that 
(£ — C)/Po corresponds to X,(0) — Xj(0) after statistical 
averaging. Numerical calculations show that this is in- 
deed valid as far as the longtime behavior is concerned, 
but transiently there arc additional modifications due 
to triple correlations such as (i/>(— p, 0)$(— q, 0)ip(—k, t)). 
Before discussing these modifications, let us develop non- 
linear closure theory for C, which gives the triple corre- 
lation as a byproduct. 



IV. NONLINEAR THEORY FOR FINITE-TIME 
EFFECTS 

A. Inclusion of nonlinearity: DIA for SFD 

The expression for 2pDC in Eq. (|3.3[) is based on the 
linear approximation to Eq. (|2.12[) , which is only asymp- 
totically valid for sufficiently long time. For finite values 
of t, there should be a correction to Eq. (|3.3|) describing 
the transient behaviors of 2pDC and MSD; to find this 
correction, now we develop a nonlinear closure theory for 
the correlation C. 

In short, what we present here and in the next sub- 
section is a systematic derivation of MCT equation for 
C. Except for some minor (but important) differences, 
this is analogous to the attempt of a number of authors 
|H E3, 5H who tried to re-derive MCT equation with the 
Martin-Siggia-Rose (MSR) formalism [49[ . They were 
obstructed by the problem of inconsistency with FDT: 
this difficulty is inevitable for dense colloidal suspensions 
or supercooled liquids, as long as one uses the standard 
MSR formalism with the plain p(k, t) as the field variable 
[5fjl - l52j | . Only some special classes of Langevin equations 
are free from this difficulty. Two such classes are known: 
one is the class of models whose nonlinearity comes from 
the gradient of the thermodynamic potential (entropy) 
alone, such as the p-spin model, referred to as "Class 
I" by Miyazaki & Reichman; and the other one (called 
"Class II" ) is exemplified by liquid models with Gaussian 
approximation to the entropy [39l |48| . Fortunately, our 
equation for tp belongs to Class I and therefore consis- 
tency with FDT is expected. 

Let us return to the nonlinear Langevin equation 
(|2.12p governing tp, with the random force statistics in 
Eq. (|2.14l) . The correlation C is then subject to an equa- 
tion containing triple correlations: 



TV 



C(k,t) 



J* E Vf ^(- P) ^(- g ,^(-M)>, 



fc+p+(j=0 



with the 0(V; 3 ) term in Eq. (f2~T2]) discarded. Note the 
absence of (/l^) on the right-hand side; this term van- 
ishes because /l is not multiplicative, which means that 
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(fL{k,t)fi,(—k',t')} is independent of the ip'a according 
to Eq. (|2.14[) . To provide a closure to this equation, we 
apply the formalism of direct-interaction approximation 
(DIA) [H-Hl [Hj]. The procedure of the calculation is 
essentially the same as that in Ref. [54| and is explained 
briefly in Appendix [B] As a result, we obtain a set of 
equations: 

( 5 * + ^ fc2 ) ^ (M) = / dt'M G {k,t-t')C(k,\t'\) 

+ / dt'M c (k,t-t')G{-k,-t') 

Jto 

(4.1) 

(d t + ^k 2 ^j G(k,t) = J dt'M G (k,t-t')G(k,t') 

(4.2) 

for the correlation C and 

C<M-0-<o<M*,0>-(|g£$). 

with to (< 0) denoting the time at which the "direct inter- 
actions" arc switched off, and M G and M G are memory 
kernels given by 

AT 2 

M G (k, s) = —J2 VfVf C(p, s)G(-q, s), (4.3) 
9T 2 

M c {k,s) = —Y,(Vn 2 C(p,s)C( q , S ). (4.4) 

Note that, in the present case, the propagator (Kraich- 
nan's response function) G is essentially equivalent to the 



response function to an externally applied probe force, 
because the random forcing term pofh in Eq. (|2.12[) is 
not multiplicative but additive. 

In regard to the choice of to , we consider two possibil- 
ities. Choosing i — > — would admit a solution of the 
form C(t) = G(t)C(0), which corresponds to the "Class 
II" approximation; we do not take this choice, as this 
would require the Langevin equation to belong to Class 
I and II simultaneously, leading to a result that is ei- 
ther trivial or inconsistent with FDT. Instead, we take 
to — > —00, so that Eqs. (|4.ip and (|4.2p become identical 
to the standard one-loop result of the MSR formalism 




B. Label-based MCT equation for SFD 

In principle, Eqs. (|4.ip and (|4.2I) with the initial con- 
ditions 

C(k,0) = ^-, G(fc,0) = l 

should suffice for determination of C and G. However, as 
soon as we start to calculate them in this straightforward 
way, we find us confronted with difficulties. Since we 
have to — > —00, the equations are acausal. Besides, the 
memory terms seem to suffer ultraviolet divergence. To 
elude these difficulties, we differentiate Eq. (|4.1[) in regard 
to t and add it to Eq. (|4.2j) multiplied by aok 2 with some 
constant a , to write an equation for d t C + a k 2 G: 



'c> 4 + ^fc 2 ) [d t C(k,t)+a k 2 G(k,t)] = I dt'M G (k,t-t')[d t C(k,t') + a k 2 G(k,t')] 



+ / dt' {M G {k,t-t')dt.C(k,-t')-[d v Mc(k,t-t')]G{-k,-t')} . 

Jt 

(4.5) 



The second term on the right-hand side includes V's through the memory kernels, which we rewrite by substituting 
Eqs. (|4.3[) and (|4.4[) . Subsequently, introducing W by V^ q — D*k 2 Wk pq and making use of the symmetry of W [see 
Eq. (|2.13|) ]. after some algebraic manipulation, we find 



[the integrand in the 2nd term on RHS of Eq. (|4.5 
~N 



ATI 

— Dlk 2 ]T W 2 pq [q 2 G(~q, t - t')d t ,C{k, -t') - k 2 G(k, -t')d t >C(q, t - t')] 



2 



^Dlk 2 Y,W 2 pq x 

{q 2 G(-q, t - t') [d t >C(k, -t') - a a k 2 G{k, -if)] - k 2 G(k, -t') [d t >C{-q, t-t')- a k 2 G(-q, t - t')] } . 

I 



Then Eq. (j4.5[) can be replaced with a simpler relation in the sense that both sides of Eq. (|4.5[) vanishes if 



d t C(k) + a k 2 G(k,t) =0 (forVfc), (4.6) Eq ' ^ holds - Thc constant is determined to be 
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«o = D*/L 2 by the initial condition. 

Taking notice of the property of Eq. (|2.12[) that the 
propagator G is equivalent to the response to the probe 
force, we note that Eq. (|4.6[) states the fluctuation- 
dissipation theorem (FDT), which can be derived directly 
from the Langevin equation (|2.12[) through the distribu- 
tion function [55l - [57j . In other words, FDT is already 
included in Eqs. (|4.1|) and ()4.2[) . This inclusion is a re- 
markable feature of Eq. (|2.12[) . or Eq. (|2.8|) . if we com- 
pare it with an analogous calculation starting from the 
Fourier representation of the "Eulerian" equation (|2.2[) . 
as opposed to the "Lagrangian" equation (|2.8j) . In the 
"Eulerian" case, the step corresponding to the rearrange- 
ment of Eq. (|4.5|) turns out to be inconsistent with FDT 
[54^ This inconsistency is due to the hidden depen- 
dence of f p (x, i) on p in its statistics in Eq. 



known 

as the multiplicative noise [39j, which makes the "Eule- 
rian" equation intractable with DIA-likc expansion. Con- 
trastively, the statistics of /l on the RHS of the "La- 
grangian" equation (|2.8|) is given by Eq. (|2.9|) which is 
independent of ip. This is why the DIA equations (|4.1[) 
and (|4.2|) successfully reproduce FDT. 

Equation (|4.6[) allows us to eliminate G from Eq. (|4.1|) 
and thereby elude the difficulties mentioned at the be- 
ginning of this subsection, as it implies 



G(k,t) 



from which we can show 



1 



Ctok 



;d t C(k,t), 



M G {k,s) 



;d s M C (k, S ). 



Then we substitute it into Eq. (|4.1j) , and the result reads 



d t + ^k 2 



1 



aok 2 



M c (k,0)C- / dt'M c {k,t-t')d v C{k,t') 
Jo 



The source of the ultraviolet divergence is now isolated in 
Mc(k, 0), which we should discard, as this term seems to 
have originated from an inappropriate treatment of the 
self- interaction in DIA J58| . Thus we arrive at the MCT 
equation: 



d t + ^jk 2 j C(k,t) = ~ J dt'M{k,t-t')d t >C(k,t') 



where 



(4.7) 



M(k, s) 



M c (fc,a) 

aok 2 
2L A 

-D*k 2 Y, W 2 qk C(j>, S )C{q,s). (4. 



N 



C. Solution to MCT equation 



Now the finite-time correction to Eq. (|3.3[) for 2pDC 
is within our reach: all we need to do is to solve 
Eq. (|4.7|) and substitute the solution C into the ex- 
tended Alexander-Pincus formula (|3.2[) . Although one 
may switch to numerical remedy, here we prefer to stick 
to the fully analytical calculation, which is possible by 
assuming the dilute limit (po<J — > +0; S = 1, D c = D). 
This does not trivialize the problem, because nonlinearity 
still exists due to 



1 



1 + V 



= 1 - v + v 



and therefore the RHS of the MCT equation (T4~71) docs 
not vanish. Let us evaluate it using the linear solution 
Eq. (|2.18|) as the zeroth approximation valid for t —> +00, 
which now reads C(k,t) ~ L~ 2 e~ D * k * as S = 1. 

To start with, we calculate M(k, s) by substituting the 
approximate solution into Eq. (|4.8[) . Parametrizing the 
variables in the summation as (j>, q) = (k/2 + m, k/2 — 
m) and denoting the wave number interval with Am = 
2tt/N, we find 



M(k,s) 



D,k 2 



^exp 



-D* 



2m' 



Am. 



The summation is then replaced with an integral, which 
readily yields 



M{k,s) 



D*k 2 
\J2i:D*s 



-.e 2 



and therefore 



[RHS of Eq. (T4T71) ] 



Dlk 4 
L 2 



-D,k A t 



df 



y/2nD*(t-t') 



-lD«fc 2 ((-f') 



(4.9) 



Though the integral in Eq. f|4.9[) can be evaluated rig- 
orously in terms of the error function with an imaginary 
argument, it is more convenient to evaluate it by expand- 
ing the integrand in powers of t' — t, as the largest contri- 
bution to the integral comes from the vicinity of t' = t. 
Thus we find 

[RHS of Eq. (TO) ] 



2f.4 



Djk 
L 2 



71- 3 V 2tt 



which allows us to integrate Eq. (|4.7[) as 
1 



l + lJh\D*tfl 2 + -? T = 

3V7T 15V27T 



k\D*tff 2 



p+q=k 



(4.10) 
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It should be possible, at least in principle, to substitute 
Eq. (gUni into Eq. flU) and the RHS of Eq. g^l for the 
second approximation, but for the present let us content 
ourselves with this first approximation to go ahead. 



using the formula (|3.2p . The procedure is analogous to 
that for the derivation of Eq. p.3[) from the linear solu- 
tion in Eq. (f2A8|) . 



D. 



Effects of the nonlinear term on transient 
behaviors of MSD and 2pDC 



Since C is now available in Eq. (|4.10|) as a result of 
nonlinear closure theory, we can evaluate (R(£, t)R(£' , t)} 



If we take into account the term of order (Z?»i) 3 / 2 and 
ignore that of order (D*t) 5/2 in Eq. (|4~T0|) . by substitut- 
ing Eq. (|4.10p into the formula (|3.2p we obtain 




cxp 



'V2 



Ap 2 Dt j 



Pi 



2^ e -6 



2p VDi 3tt 



V2 _ 2 
Po 



l\2- 



2p 2 Dt 



A2" 



exp 



4/^ 



(4.11) 



where = (£—£') / {2pov Dt) and the function ip is defined 
in Eq. (|3.5p with 5 = 1. As a special case for £ = 
Eq. (|4.1ip gives correction to (i? 2 ) oc Vt: 



/ P 2\ 2 [Dt V2 _ 2 

PO V 7T 37T 



(4.12) 



The first term reproduces Eq. (|1.7p with S = 1, while 
the second term gives a correction to it. The contribution 
from the higher-order terms in Eq. f|4. 10[) slightly enlarges 
the coefficient of the correction term, but the form of 
Eq. (|4.12p itself is not affected. 

Without the second term in Eq. (|4.12l) . plotting (R 2 ) 
against s/t would yield a graph of a straight line pass- 
ing through the origin. In actuality, the second term in 
Eq. (|4.12p shifts the asymptotic straight line, making a 
positive intercept on the yi-axis and a negative intercept 
on the (i? 2 )-axis. These intercepts are already visible in 
the inset of Fig. [IJ and also in Fig. [5] (a plot analogous to 
Fig. [T]but with the axes rescaled, as we explain below). 
The presence of the correction itself is probably not sur- 
prising, because the short-time diffusion should behave 
as (R 2 ) oc Dt before the collisions begin to take effect; it 
is more noteworthy that, since this single-particle behav- 
ior plays the role of the "mode-coupling" in the Fourier 
representation, the description of the transient behavior 
requires a nonlinear theory such as MCT. 

The prediction of the nonlinear theory in Eq. (|4.12l) . 
including its dependence on the density po, is compared 
with the computed MSD in Fig. [5] As the quantitative 
comparison requires us to take into account the effects 
of the finite density, we revived S = 5(0) in the first 
term of Eq. f|4.12[) . and plotted p\ (R 2 ) against Sp \JD c t. 
The MSD computed for three different values of density 
(p a a = 1/4, 1/8, and 1/16; see the inset) seem to collapse 
into a single curve whose asymptote is the straight line 




1 2 3 



FIG. 5: (Color online) Comparison of Eq. (|4.12[) with numeri- 
cal data, by means of rescaled plotting of MSD versus time for 
three different values of the density: po = N/L = (l/4)cr _1 , 
(1/8)<t -1 , and (1/16)<t _1 . The number of the particles is fixed 
at N = 2 15 = 32768. The solid line represents the prediction 
of the nonlinear theory in Eq. (|4.12|l , which is compared with 
Kollmann's law (|1.7[1 in the dotted line. The inset shows the 
same data without rescaling, using the simple nondimension- 
alization with a 2 and a 2 / D. 



given by Eq. (|4. 1 2[) . Improvement of Eq. (|4. 1 2[) account- 
ing for the small deviation from the straight line will be 
performable with a careful numerical calculation of the 
MCT equation, which will be reported elsewhere. 
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b 




X,(0)]/A(i) 



FIG. 6: (Color online) Eq. (|3.3[) compared with simulation 
data (with the same N and po as in Figs. [T] and [4}, on the ba- 
sis of the initial distance Xj (0) — Xi (0) instead of the particle 
numbering. The solid (red) line shows the theoretical predic- 
tion for t — > -(-oo by Eq. (|3.3|l . while the simulation data are 
plotted with the same symbols as in Fig. [4] 



V. OTHER FORMS OF FOUR-POINT 
CORRELATION DERIVED FROM 2PDC 



Behavior of 2pDC as a function of the initial 
distance 



As we have already shown in Fig. [H the theoretical 
prediction for (R(£,t)R(£' ,t)) is almost perfectly consis- 
tent with numerical calculation. However, taking the 
label distance £ - £' or 9 = (£ - £')/\po\(t)] for the 
horizontal axis of the graph is sometimes inconvenient, 
for example, when we try to compare the result with 
two-dimensional or three-dimensional numerical calcula- 
tions. For this purpose, it may be more convenient to 
re-express the result as a function of the initial distance, 
Xj(0) — Xi(0), and plot the two-particle displacement 
correlation against [Aj(0) — Xj(0)]/A(i). Such a graph is 
shown in Fig. [HJ The analytic curve in Fig. [5] is drawn by 
simply equating [Xj(0) -X;(0)]/A(t) with 6 in Eq. (|3~5]) . 
This seems to be valid for t — > +oo, but a considerable 
discrepancy is seen for shorter times. 

Although Xj(0) — Xi(0) and (Sj — Sj)/po are equal 
on the average, generally they are different, as is evident 
from the relation 



X J (0)-X i (0) = 



- 1 + V>(£,0) 
Po 



d£. 



(5.1) 



This difference is responsible for the discrepancy in Fig. [6] 
for finite t. Taking this difference into account, we can 
evaluate 2pDC theoretically as a function of the initial 
distance. Although its complete evaluation is out of 
the scope of the present paper, as it seems to involve 
wavenumber integrals that are difficult to perform ana- 



lytically, we can outline the procedure of the calculation 
at least. 

With the value of the initial distance denoted with 
d, the function that gives 2pDC, which we denote with 
Xft(d,t), is formally written as 

X R (d,t) 

= \X ff 5 ( Xl ~ X2 -d)fl(£i,*)fl(6,*)dxidEa 



(5.2) 



where 

Xi = Xi(£i) = x(£i,0), x 2 = £2(6) = x(£ 2 ,0), 
and therefore 



x 2 - Xl 



1 

Po 



6-6+/ ^(£,o)d£ 
6 



Changing the variables of the integral in Eq. (|5.2p from 
(xi,x 2 ) to (£1,62) witn 

1+^0) 1 + ^(6,0) ^ 

dxi = d£i, dx 2 = d£ 2 , 



Po 



and introducing the Fourier representation of the delta 
function, 



6(x) = jY / e- iq * 



( 1 



L^" \2tt/L 
we rewrite Eq. (|5.2p as 

where k = q/ po and 



(5.3) 



R(k,t) = - e iqx R(£(x,0),t)dx 
L Jo 



NJ 



N 



cxp ifc 



£ + / Vd£ 





i?(£,i)(l + V)d£, 
(5.4) 



with 



^ = V(£,0), / V>d£ = / ^(£,0)d£. 
Jo Jo 



Then we express ip i n Eq. ()5.4[) with ip in Eq. (|2.11| and 
also substitute Eq. (|3.1j) . After some rearrangement, we 
obtain 

- l g(M)-j(M) 

X =ijfe 

X E ^,0)[^)-^-,0)] +o( -.3^ (55) 



A 



ip 
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Substituting Eq. (|5.5[) into Eq. (|5.3j) yields an expression 
of Xr(cI, t) that consists of two parts: the first part sim- 
ply reproduces Eq. (|3.2[) with £ — £' replaced with po^, 
and the second part involves triple correlations such as 
(^(-p,O)i/i(-g,i)0(-M)) with p + q- + ft = 0. These 
triple correlations can be calculated with DIA, and as a 
result, we obtain a correction term whose magnitude rel- 
ative to the leading term decreases in proportion to t~ x l 2 
for t — > +00. Detailed results of the calculation will be 
reported elsewhere. 

We note that the definition of Xr in Eq. ()5.2[) is readily 
generalized to three-dimensional cases, as 



X(d,t) = (jp / 



L 3 



( TV 2 



» j 



52 (ry) 



Ri ® R 7 



(5.6) 



where ry = r,- - r 4 , Rj = R(&, t) = r^t) - r*(0), and 



i' i' 



This X is similar to the quantity calculated by Donati e< 
al. [37[ | ( <7 U in their notation) , except for two main differ- 
ences: in their g u , the two-body density (72 is absent in 
the denominator, and a product of scalar displacements, 
|Rj| |Rj|, is used instead of Ri <E> Rj. The presence or 
absence of gi is not essential, though it indeed makes 
it difficult to define X for small values of the initial dis- 
tance, |d| < a, in which we are not interested. The other 
difference is crucial: the absolute value signs obstructs 
analytical evaluation of g u even in the one-dimensional 
cases. Besides, the tensorial character of X can provide 
a useful information on the geometry of the collective 
motion in the three-dimensional glassy systems. We will 
return to this point in Sec. IVI1 but before that, let us 
relate the one-dimensional 2pDC to Xi- 



B. Calculation of X4 from 2pDC 

With the knowledge of the displacement correlation 
(R(£, t)R(£',i)) in Eq. (|4.11[) . we can also calculate a one- 
dimensional version of a quantity which is commonly re- 
ferred to as X4 W- To be precise, we consider the Q-based 
X4 [IH EH j as opposed to other variants of %4 such as the 
.F-based Xi defined through the fluctuation of the 
intermediate scatterin g fu nction F or its self part. If we 
consult Glotzer et al. [13| and adapt their equations (4) 
and (5) for one-dimensional cases, we have 



* 3 



l (Q 2 )-(QY 

kv.T N 2 



(5.7) 



(5.8) 



This type of four-point correlation function has been 
studied by many authors [HI Q f60M63| as an indicator 
of cooperative motion in glassy systems. To our knowl- 
edge, most of these studies are based on direct numerical 
simulations of particle systems and there are also exper- 
iments grounded on observation of particles, but analyt- 
ical calculations are quite rare. What makes it difficult 
to calculate Xi analytically is that, in the usual formula- 
tion, the four-point correlation implies four-body corre- 
lation, i.e. it involves four bodies. More concretely, as Q 
in Eq. (|5.7p already contains double summation, calcu- 
lation of X4 requires dealing with quadruple summation 
whose summand involves four particles simultaneously; 
this would be a hopeless task. 

To facilitate calculation of four-point correlation, here 
we introduce two modifications to Eqs. (|5.7[) and (|5.8|) . 
Firstly, we target on the "self part" (i = j) of Q and its 
contribution to X4j denoting them as [olj 



L (Ql)-(Qs) 

k n T N 2 



(5.9) 

Since Glotzer et al. [l3| reported that the contribution of 
the self part (i = j) is dominant over that of the distinct 
part (i 7^ j), it is justifiable to calculate xf instead of %4- 
Secondly, as the overlap function S a , we adopt a Gaussian 
function 



2 i 2 

-r j a 



(5.10) 



instead of the step function used by Glotzer et al. fl3j |. 
We note that, although there exists a variant of %4 from 
whose definition the probe length can be totally expelled 
[46| . the probe length a is indispensable to Qs- 

From Eq. (|4.1ip we already know the covariance 
(RiRj) for all and for arbitrary t (within a certain 
limitation, of course). The problem is how to evaluate 
(Qs) an d (Qs) m Eq- (|5.9[) using this information. This 
is possible, if we assume that (i?i, R2, . . . , Rn) is subject 
to a joint (multivariate) Gaussian distribution, which is 
determined uniquely as the covariance is given and the 
mean is known to vanish. For the purpose of calculating 
(Q|); ^ suffices to determine the two-body distribution 
function for (i?i, Rj), which we denote with 



P(Ri,Rj) 
1 



: exp 



(R 2 ) (R 2 + R 2 ) - 2 (R.Rj) R^ 
2A„- 



wherc 



A, tJ = (R 2 ) 2 - (R lRj )\ (R 2 ) = (R 2 ) = (R 2 ) . 

Using this joint distribution function P(Ri,Rj) and 
adopting Eq. (|5.10[) for the overlapping function, we ob- 
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tain 



(Sa(Ri)) 



and 

(S a (Ri)S a (Rj)) 



~5 a {R l )P{R l ,R j )dR l dR J 
1 



(5.11) 



SaiR^SaiR^PiR^R^dRidRj 
1 



(5.12) 



note that Eq. (|5. 12|) is confirmed separately for i ^ j and 
i = j. With Eq. (|5.1ip and (|5.12|) . now we can evaluate 
X4 in Eq. (|5.9j) . taking the uniformity of the system into 
account. As a result, we obtain 



xl 



N 2 k B T 



* 3 



L 



Nk B T 



E 



\ 



2<i? 2 )\ 4(i? 4J R w ) 2 1 + 



2<i? 2 ) 



(5.13) 



Note that the double summation J2j (Sa(Ri)S a (Rj)) 
in Eq. (|5.13[) is a result of the simplification by the re- 
placement of X4 with xt (retaining only the self part): 
if this simplification were not introduced, we would have 
to struggle with a quadruple summation such as 



whose evaluation would be much less workable than 

(SaiRMRj)). 

Before applying Eq. (|5.13|) to SFD, we can test it with 
free Brownian particles. From the Langevin equation ob- 
tained by setting V = in Eq. (|1.5p . we have 



(R 2 ) = 2D[t-T B (l-e-t/™)] (i=j) 




{RiRj) to a^j) 

(5.14) 

where t b = m//j,. This is substituted into Eq. (|5.13p . 
which yields 



X4 S 



PoksT \^1 + 4(R 2 ) /a 2 

(xl)solo 



1 + 2 (R 2 ) /a 2 



(5.15) 



for free Brownian particles; note that all the contribu- 
tion comes from the term with I = in Eq. (|5.13|) , which 
we refer to as the "solo" part. Taking notice of the t- 
dependence of (R 2 ) in Eq. (|5 . 1 5[) and making some cal- 
culation, we find (xi)soio to have a peak at the instant 
when (-R 2 ) = 2.6 a 2 approximately; see the short-time 
side of Fig. [TJa). Obviously, this short-time peak is ir- 
relevant to particle interaction. After this peak, (xl^soio 
decreases monotonically toward zero, in proportion to 
£-1/2 f or t +ocx 



Now let us calculate xl for SFD, combining Eq. (|5 . 1 3[) 
with the result for (RR) in Eq. (|4.1ip . We evaluate xf 
in Eq. (|5.13[) as a sum of the "solo" part (I = 0) and the 
collective part (contribution from the terms with I =/= 0). 
The solo part is given by Eq. (|5 . 15[) and depends on 
(R 2 ) alone, for which we use Eq. (|4.12p that was ob- 
tained by setting £ = £' in Eq. (|4.1ip . As a matter of 
course, we must exclude cases of very short time, for 
which Eq. (|4. 12[) predicts (R 2 ) to be negative; this is out 
of the validity range of Eq. (|4.12p . Subsequently, to eval- 
uate the contribution from the terms with I ^ 0, we use 
the asymptotic form of Eq. (|4.11l) for t +oo, expressed 
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(b) 0.6 




r l/4 

FIG. 7: (Color online) Short-time and longtime behavior of xl for SFD, compared with the analytical prediction given by 
the sum of (x4)soio in Eq. (|5.15[) and (xl)coii m Eq. (|5,17|l . with the effect of S ^ 1 and D c /D 7^ 1 taken into account. The 
probe length (radius of the overlapping function) was chosen as a = 0.5 a. The time is non-dimensionalized with a 2 /D. (a) 
Short-time behavior for the case with N = 256 and po — N/L = (1/4) a -1 . The computed data are plotted with circles and 
the analytical prediction is shown with a solid (red) line, (b) Long-time behavior for po = (1/4)<j -1 (plotted with circles) and 
for po = (l/16)cr _1 (plotted with crosses). The solid lines show analytical prediction for the two cases. Note the asymptotic 
behavior of the curves that takes the form of a straight line in the graph, exhibiting the decay of (xl)soio in proportion to i~ 1//4 
and the finite value of (xf)coii that remains for t — > +00. 



as a self-similar solution in Eq. (13. 5]) : 



by rewriting it as 



(RiR i+l ) * 0, = — ^ 



Po 



PoA(i) 



The collective part is thereby written as 



(Xl)coll 



I y 

>/r2\ 



p k B T 1 | 2<g) 1M 
a 2 



<p(0) 



The number of the particle contributing the the sum is 
estimated to be N co i\ ~ 1/A9 = poX(t), which yields, 
approximately, 



(xl) 



1 



coll 



Po X(t) 



p k B T 1 + 2 (R 2 ) /a 2 
a 2 X(t) a 2 p 
k^T X Jr 2 ) ~ feT 



(5.16a) 
(5.16b) 



for £ — > +00. More precisely, the sum can be evaluated 



5>0 = 2poA(i)X; 

Z^O i=l 



2p A(i) 



1 



1 - 



V(0) 
1 



a/i-tt^)] 2 
and using the numerical value of the integral 



\ A9 



1 V d0 



/o i ^1 - 7T [¥>(£>)] 

thus we obtain 

1 poA(t) 



1 > d6» = 0.364124; 



(X/Jcoll 



x 0.364124 



p k B T (R 2 ) /a 2 

— — x V^pW x 0.364124 
pok B T 



= 0.6454 x 



k B T 



(5.17) 



and xl = (xD solo + (xf)coii, with ( xpsoio given by sub- 
stituting Eq. dUTH into Eq. (|5TT5]) . If we take S into 
account according to Eq. (|3.5p . the right-hand side of 
Eq. (f5TT7|) is multiplied by S _1 . 

The analytically calculated xf and its numerical values 
are compared in Fig.[7J The effect of finite po that makes 
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S and D c /D different from unity is taken into account 
(see Table U in Appendix |A"]) . The peak in the short-time 
regime has nothing to do with the slow dynamics, as it 
appears even for free Brownian particles. After this peak, 
(xl)soio decreases slowly, asymptotically in proportion to 
i -1 / 4 , while, reflecting the endless growth of the cluster 
size, (xf )coii remains finite for t — > +00. The behavior of 
the numerical solution is consistent with this analytical 
prediction. 

The limiting value of xt f° r t —> +00, given by (xl)coll 
in Eq. (|5.17[) , is an increasing function of the density 
po . This conclusion remains unchanged also if the effect 
of S =/= 1 is included, because 1/5 is also an increasing 
function of p . 



VI. DISCUSSION AND CONCLUDING 
REMARKS 

We have studied one-dimensional systems of Brown- 
ian particles with repulsive interaction, regarding it as 
a simplified model of the cage effect. On one hand, 
the cage confines every particle in a narrow space. On 
the other hand, since the particles are mutually caged 
and therefore forbidden to move uncoopcratively, they 
must either wait still or move together. Thus the cage 
effect involves correlated motion of numerous particles, 
which is visualized as dynamical clusters in space-time 
diagrams (Fig. [3]), with the diffusive correlation length 
A(i) = 2\/D c t. The substance that diffuses is not the 
particles but the space between them, represented by tp 
in our theory. The two-particle displacement correlation 
(2pDC) is shown to be a useful indicator of the correlated 
motion. Asymptotically, 2pDC becomes self-similar: it 
suggests a matryoshka-like structure, such that the small 
cages are confined in larger cages, which, in turn, are 
caught in still larger and slower cages. 

Using the analytical result for 2pDC which is valid 
both transiently and asymptotically, we demonstrated 
how to calculate the Qs-based %4 (denoted with xt)- 
Despite the endless growth of X(t) and the absence of 
a relaxation, the result in Eq. (|5.17j) shows that xt con ~ 
verges to some constant for t — > +00. An implication of 
Eq. (|5.17[) is that xf , and probably %4 in general, does 
not give a straightforward representation of the cluster 
size. Indeed, X(t) is in the numerator of Eq. (|5.16p or 
(I5.17|) . but the result is modified by the denominator, 
or a prefactor 1/(1 + 2 (R 2 ) /a 2 ) originating from (Qs) 2 , 
which cancels the temporal growth of the cluster size. 
In three-dimensional systems, a direct relation between 
X4 and the cluster size is expected only for some limited 
time scales corresponding to the plateau of the MSD. 

The adoption of the Lagrangian description enabled us 
to reproduce Kollmann's asymptotic law for MSD within 
the liquid-theoretical framework, calculate 2pDC analyt- 
ically, and develop the "Lagrangian MCT" without vio- 
lating FDT. The Lagrangian MCT equation was solved 
approximately by analytical means, giving a finite-time 



correction to the asymptotic law. We have implemented 
the Lagrangian description by explicitly introducing the 
label variable £ and thereby constructing a stretchable 
coordinate system that sticks to the cages everywhere. 
Probably some aspects of glassy dynamics, such as dy- 
namical heterogeneity characterizable by bond breaking 
[H, HH, [65j], may require the Lagrangian description by 
nature when its continuum counterpart is sought . 

The Lagrangian description in higher dimensions may 
not be so simple as in one-dimensional cases, but it is not 
impossible. In three-dimensional cases, a triplet of label 
variables (£, 77, £) is expected to be related with p and 
Q by equations analogous to Eq. (|2.5|) ; see Eqs. (6.6) 
and (6.7) in Ref. [Hj]. Besides, we could adopt some 
methods from three-dimensional theories of turbulence in 
which Lagrangian correlations are used [30l - l32j . Turbu- 
lence theoreticians have even considered the Lag rangian 
dynamics of a tetrad (four material points) [661 ] . whose 
two-time correlation involves eight points in the space- 
time. 

In contrast to the "Eulcrian" (standard) MCT in which 
cage effects are represented by the memory kernel (not 
successful in SFD), the Lagrangian theory can dispense 
with the memory integral as far as the asymptotic behav- 
ior is concerned. A pivotal role is played by the modified 
and extended Alexander-Pincus formulae in Eqs. (|2.15[) 
and (|3.2[) . Its linear version, namely Eq. (|2.17p or its mul- 
tidimensional extension, has been used in the context of 
glassy dynamics by several authors l43j , who limited 
themselves to the approximation with linear elasticity. 
Since the Eulerian and the Lagrangian variables are ap- 
proximately interchangeable in the description of small 
elastic deformation, they did not bother to distinguish 
the two descriptions. Needless to say, this treatment fails 
when the system is more liquid-like. In an attempt to in- 
troduce the a relaxation into the calculation of xa based 
on the "elastic" theory, Toninelli et at feared that it 
would make the model inconsistent, because the underly- 
ing lattice, needed to define the deformation field, would 
be totally melted Probably this is too pessimistic: 
the "melting" of the lattice does not make the theory 
totally inconsistent but requires more careful distinction 
between the Eulerian and the Lagrangian coordinates. 

There will be another modification to the theory of 
Toninelli et al., when departing from the linear elasticity 
and trying to consider liquid-like behavior. Their formula 
corresponding to Eq. (|3.2[) . namely Eq. (A3) in Ref. [f|, 
reads 

/ \ r 1 - p- Dk2 ' 

(i?(d)i?(0)) oc / - 2 e- ik d d' ld k (6.1) 

in our notation, as they seem to have identified the 
Langevin equation for the displacement field with the 
Edwards- Wilkinson equation [42( • In regard to Eq. (|6.ip , 
we suspect that the vectorial character of the displace- 
ment is not adequately taken into account. Probably one 
needs to decompose R into the longitudinal and trans- 
verse components, and tract them more carefully. Unlike 
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the two sound modes in elastic solids, the two modes 
in the liquids can have quite different nature: the liq- 
uid may resist compression strongly but the resistance to 
shear may be much weaker. 

In Eq. (|5.6|) , we have proposed to define the three- 
dimensional 2pDC as a tensorial quantity X. Due to the 
isotropy of the system, X must be a sum of the longitu- 
dinal and the transverse components: 



X = X 



d(g)d 
d 2 



1 - 



d® d 

d 2 



It is quite likely that X« and X± will be character- 
ized by different correlation lengths. Taking the two 
different correlation lengths into account, we can ex- 
tend the present theory phcnomcnologically to the three- 
dimensional cases. From the inferred distribution func- 
tion for the displacements of two particles, shown in Ap 
pendix [CJ we can calculate 
Subsec. EE as 



xf in the same way as in 



1 



(l-«) 2 A|,Ai 



knT 



1 



2A° 



(6.2) 



where A|| = A||(i) and Aj^ = X±(t) denote the two cor- 
relation lengths, X° = X°(t) is related to the MSD of 
the caged particles, and a = a(t) stands for the relative 
number of the particles that have hopped. Assuming the 
i-dependence of these four quantities phenomenologically 



as 



<Dt 



A± = (l + v 7 ^) 4, 



An = 



Dt 



l + Dt/£ 2 ' 



, t 

a = tanh — , 



we can plot xl against t for some different values of the 
time scale r. The curves in Fig. [5] reproduce some basic 
features of the Q-bascd X4 calculated by Lacevic et al. 
[l~H . such as the shape of the uphill that looks steeper 
near the peak in this semi- log plot. 

To go beyond the linear theory, nd-dimensional ver- 
sions of the extended Alexander-Pincus formula should 
be developed. We suppose that the nd-dimensional for- 
mula will contain, instead of C, correlations of deforma- 
tion tensor. This is not so formidable as it may appear, 
because many components of the correlation tensor will 
turn out to vanish or to have the same value as some other 
component. The formula is now under development and 
will be reported elsewhere. 

It is also noteworthy that, in the derivation of MCT, 
the difficulty of FDT violation disappeared quite natu- 
rally due to the Lagrangian description. As is pointed 
out by Miyazaki & Reichman [39j . it has been difficult 
to construct realistic models which does not violate FDT 
and can incorporate the effect of structural changes em- 
bodied in S(k) at the same time. Since the "Lagrangian 



x 




Dt/tl 



FIG. 8: (Color online) Behavior of \4 given by Eq. (|6.2|) . 
calculated through the three-dimensional two-body displace- 
ment distribution function by phenomenologically extending 
the present theory to the three-dimensional cases. The three 
curves correspond to different values of r: from left to right, 
t = W 3 el/D, 10 4 el/D, and lO^o/D. 



MCT" is now shown to be consistent with FDT, study of 
its behavior for different S(k) may be quite intriguing. 

For possible extensions in the future, we can mention 
several directions. For example, one may include weak 
attractive interaction and analyze the effect of the change 
in S(k) on the transient behavior of SFD. One may also 
study nonequilibrium behavior by driving the particles 
with an external force or changing the temperature sud- 
denly. Another interesting proposal is to permit over- 
taking as a rare event, which may play the role of the a 
relaxation. We performed simulations with some finite 
interaction potential, and a preliminary result juTj shows 
that normal diffusion is observed for y max = ksT, while 
for V max = 5fceT the behavior is essentially that of SFD 
in the time scale of the simulation. The problem is to 
make a theory that can handle the crossover between the 
two limiting cases. The theory allowing for rare overtak- 
ing events may bridge the gap between the purely one- 
dimensional SFD and the behavior of three-dimensional 
rod polymers [TH, [68| . 

We could also study double-file diffusion, which would 
be analogous to two-lane models of traffic flows. If a "lane 
interaction" is also introduced, the system would have 
also something common with the Matsukawa-Fukuyama 
model of friction formulated on a ladder lattice [63, [7(| . 
The study of the double-file diffusion may shed light to 
many related systems in which frustration is dynamically 
created and annihilated, such as template-assisted pat- 
tern formation of colloid particles on a substrate with 
parallel channels [7l|, frustrated Josephson-junction ar- 
rays in a magnetic field [lO, [72[, and — hopefuly — also 
three-dimensional dense colloidal suspensions. 

T he p resent work, in combination with the previous 
one (24[, is intended as several first steps toward a fu- 
ture theory of three-dimensional glassy liquids, which 
will make it possible, for example, to replace the semi- 
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phenomenological curves for %| in Fig. [8] with a first- 
principle theoretical calculation. Although the present 
theory is still embryonic, it already suggests that one 
of the important ingredients of the future theory may 
be the displacement distribution function of two or more 
particles. In the case of computational analysis, probably 
we should not insist on some favorite statistical quantity 
alone, nor content ourselves with the single- particle van 
Hove function, but try to deduce some suitable distribu- 
tion function behind the computed statistical quantities. 
The present analysis of one-dimensional cage dynamics 
and the concepts used for it will provide a useful frame- 
work both for analysis of numerical data and for new 
development of theory of glassy liquids. 
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Appendix A: Numerical calculations 

Here we describe how we integrated the one- 
dimensional Langevin equation (jl.5[) . and also how we 
evaluated the statistical quantities, including (R 2 ), D c , 
and S, from the numerical solution. As the system con- 
tains N particles in a periodic box of the size L, the mean 
density is given by po = N/L. 

The potential V in Eq. (|1.5|) was specified as 

V(r) = l V ™{ 1 -v) 2 (Al) 
lo (M>ff) 

with Vmax S> fceT. In the present numerical calcula- 
tions we adopted the value y max = 50 fceT, which is 
high enough to forbid the overtaking of the particles com- 
pletely. 

The random forcing is the zero-mean Gaussian noise 
whose variance is given by the one-dimensional version 
of Eq. (TO]), 

(f i (t)f j (t'))=2Dd ij 6{t-t'). (A2) 

Computationally, the delta function in Eq. (|A2[) was dis- 
cretized with the time interval Atf, and as the values of 
(/o, /i, . . . , /jv-i) for each time interval, N independent 
Gaussian random numbers with the variance 2D /Atf 
were generated with the Mersenne twister and the Box- 
Mullcr transform. 

With V and /, given as above, in nondimensionaliza- 
tion of the governing equation using a, a 2 / D, and m 



as the units of length, time, and mass, there appears 
a nondimcnsional parameter specifiable as the ratio of 
tb = m/fM to the time unit a 2 /D, which we chose as 
1:1. Then the time integration of Eq. (|1.5[) was per- 
formed with a Verlet-like scheme. The time step At was 
taken equal to a hundredth of the time unit a 2 /D, and 
the renewal interval of the random forcing, Atf, was cho- 
sen to be the same as the time step: 

2 

At = At { = 10~ 2 x (A3) 

We also tested some different choices of At and Atf, such 
as 

^^=(10- 3 ,10- 3 ) or (10- 3 ,HT 2 ). 

To bring the system into equilibrium, we started each 
calculation at t = —T w , introducing a sufficiently long 
waiting time T w (typically T w = 10 a 2 /D), and waited 
till t = 0. Then, from the simulation data recorded for 
< t < tmaxj we calculated a desired statistical quantity 
as the average for n samples extracted from the data by 
a time shift. For example, (R 2 ) is calculated as 

1 n-lJV-l 

(R(tr-) = ^EE + o - x mf , (A4) 

i=o i=a 

where ti denotes the starting time of the Z-th sample. 
Care must be taken so that the maximal value of t' in 
Eq. (|A4[) . which equals the span of each sample, should 
not exceed the waiting time 7^; a result for a longer 
span will expose insufficiency of equilibration. Typically 
we chose t max = 5T Wl n = 100, and t\ = Z(t max - T w )/n, 
allowing the samples to overlap. 

The computation of Xn{d, t), shown in Fig.[6j was per- 
formed with discretization of Eq. (|5.2[) in which the delta 
function was approximated by a statistical bin 0.5 er in 
width. After recording Xi(t{) for every particle, we clas- 
sified every pair into a statistical bin according to 
the "initial" distance Xj(ti) — so that the K-th 

bin contains the pairs for which \Xj(ti) — X,-(t() — d K \ is 
smaller than the half width of the bin. Then we calcu- 
lated Xn(d K ,t) as the average of Ri(t)Rj(t) for the K-th 
bin, where Ri(t) = Xi(ti+t) — Xi(ti). If, instead, the sum 
of the absolute values of the data in each bin were cal- 
culated, this would be analogous to the quantity studied 
by Donati et al. [37| . 

The collective diffusion coefficient D c is determined 
from the temporal decay of the dynamical structure fac- 
tor F(q,t) [26}. We computed F(q,t) for < q <C p 
and made a linear fit for log F(q, t) = log S(q) — D c q 2 t to 
obtain the values of S(q) and D c as the fitting parame- 
ters. After taking an average for several small values of 
q, the results are summarized in TablelU These values are 
used in evaluation of the theoretical predictions, such as 
Eq. (|3.3[) . and also for rescaling of the horizontal axis in 

Fig.m 
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TABLE I: Numerical values of .5(0) and D c /D computed for 
three different values of the density. 



Po 




(i/sy- 1 




1 — 2p <7 


0.500 


0.750 


0.875 


5(0) 


0.624 


0.787 


0.888 


D c /D 


1.59 


1.27 


1.12 


S(0)^/D c /D 


0.79 


0.89 


0.94 



Appendix B: Direct-interaction approximation 

The main idea of DIA [29|, [HH for evaluation of the 
triple correlation 

(i>(-p, m-q, t)$(-k, 0)) = (j>(p, t)i>(q, t)i,(k, 0))* , 

with the asterisk denoting the complex conjugate, is to 
utilize the property of V = (V^ 7 ) that, out of its N 3 
components, almost all are zero because the condition 
Q + /3 + 7 = 0is not satisfied. The nonzero components 
of V constitute what we call triad interactions: rewriting 
Eq. (|2TT2"f as 

{d t + ^ a )i>(a, t) = ]T KH* 08. *W (7, t) + poh (a, t) 

makes it clearer that V^ 7 engages in connecting the 
"triad" that consists of a, /3, and 7. If we visualize each 
triad interaction as a triangle on a graph, the property 
of V is such that no triangle shares its side with other 
triangles. Therefore, if we "switch off" a single triad, 
say, {p,q,k} (with which we mean Vf 9 , Vf , Vf, V*« , 
Vg P , and V|*; note the symmetry in regard to the inter- 
change of the superscripts, V^ q = V| p etc.), the direct 
interactions between the three modes p, q, k are lost. 

To concretize this idea, let us suppose that an artificial 
forcing term 

(2v p k q r( P ,t)r(q,t) ( a = k) 

I = - 6 (t-t)—x }W* k r(q,t)r(k,t) (a=p) 
[ 'L2 }2Vpr(k,t)r(p,t) (a = q) 

[ (otherwise) , 

designed to cancel a single triad {p,q,k}, is applied to 
the system given by Eq. (|2.12ip . We denote the solution 
to this artificial system with -0 O = ifio({p, q, k}; a, t). One 
of the two main assumptions of DIA is that the three 
selected modes, namely ifto(p,t), ipo(q,t), and ipo(k,t) in 
this case, become uncorrelated, since the forcing I cancels 
the direct interactions. On the other hand, I is regarded 
as a small perturbation, because it cancels only a single 
triad interaction and there remain still a large number of 
triads connecting, say, p and q indirectly. Therefore the 
difference ipi = — ifio is assumed to be small, which is 
the second main assumption of DIA. 

Due to these assumptions, the triple correlation is ex- 



panded as 

(^(p,t)$(g,t)$(M)) - (&(p, *)&(M)) 

+ (&(p,f)$>(g,t)&(M)) 
+ (Mp,t)Mq,t)Mk,o)) 
+ (Mp,t)Mq,t)Mk,o)) 

+ 0($), (Bl) 

and the zero-th term, (i[>o(p, t)ipo(q, t)il>o(k, 0)), vanishes. 
Since 7 is a small perturbation and is a response to it, 
formally ipi can be expressed in terms of the propagator 
G as 

ipi(a,t) = - [ di'VG(a,<; a', t')I a > (*') (B2) 

J to 

for t > to. We substitute Eq. (|B2[) into each term on the 
right-hand side of Eq. (|B1[) to find, to our surprise, that 
the result is naturally factorized due to the assumption of 
DIA that ipo{p,t), ipo(q,t), and ipo(k,t) are uncorrelated. 

By applying the decomposition <j> = ipg + -0i to each 
triple correlation term in the equation for dtG, we are 
led to Eq. (|4.ip . Similarly, the equation for dtG contains 
(tpoGY which is evaluated with the aid of the DIA de- 
composition of G, resulting in Eq. (|4.2[) . For more details, 
see Refs. (53l.l54|. 



Appendix C: Three-dimensional calculation of \t 

Here we outline how to calculate \t by phenomeno- 
logically extending the present theory to the three- 
dimensional cases. The displacement correlation tensor 
X is related to the distribution function P(d; R i; Rj) by 

X(d) = Jj ' R,® R i P(d; R;, R J )d 3 R i d 3 R 3 ; 

in what follows, the d-dependence of P is taken for 
granted and therefore omitted. The relation can be in- 
verted if the functional form of P is known. In particular, 
if a multivariate Gaussian distribution (which we denote 
with Po) is assumed, P = Pq can be factorized as 

P (R h Rj) = P,| (r| , r]|)p±(r^ R|) (CI) 

by splitting R into the longitudinal and transverse com- 
ponents as R = R 11 + R ± (so that RH || d and R^ _L d). 
Then we introduce X° such that 

X\\(d = 0) = X±(d = 0) = X° 
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and write the two factors explicitly as 



P\ = 



2tt^ 
1 

(27r) 2 Aj 



cxp 



exp 



2A,| 



2X±(d)'Rf- ■ R^ 



2A ± 



where A,, = (X°) - X {[ (d) 2 etc. 

Then the calculation of xt will be carried out in a way 
analogous to the ID cases. Subsequently, assuming that 
two correlation lengths, Am = X\\(t) and Aj_ = Aj_(t), can 
be introduced so that 



X 



X°(i)$||(J/A||(i)), 



we estimate the number of the particles contributing to 
the sum as A^ co n ~ p^X^X 2 ^, which leads to 

1 



x2 



A n A i 



1 



2X l 



as a three-dimensional counterpart of Eq. (|5.16al) . 

Of course, Eq. (|C2[) needs to be modified by taking a 
relaxation into account. We introduce a = a(t) denoting 



the relative number of the particles that have hopped by 
the time t, and assume the distribution function in the 
form 

P(Ri, Rj) 

= Ul - a) 2 P + 2o(l - a)P x + a 2 P 2 (i ^ j) 
\(l-a)P Q +aP 2 (i = j), 

(C3) 

where Pq, governing the pairs of caged particles, is given 
by Eq. (jCljl . Note that, by integrating P(Rj,Rj) in 
regard to the second argument Rj, Eq. (|C3j) reduces to 
the van Hove function in the form 

P(R) = (1 - a)Pcagc(R) + aPhop(R). 

This implies that a = a(t) can be determined, in prin- 
ciple, as a fitting parameter for the van Hove function. 
If we assume, for simplicity, that the correlation of dis- 
placements is totally lost after the hopping, we have 



(C2) Pi(R,R,) 



Plop (R? )Pcagc (Rj ) ~l~ Pcagc (Ri )Phop (Rj ) 



P 2 (Rj,Rj) — Phop(Ri)Piop (Re- 
using the distribution function in Eq. (|C3[) supplemented 
with the above expressions, we evaluate xf an d obtain 
Eq. (IQj) . 
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